Method and apparatus for assessment of fluid responsiveness

ABSTRACT

Methods and apparatus for determining a cardiac parameter from cardiovascular pressure signals including arterial blood pressure (ABP) and the photoplethysmographic signal to quantify the degree of amplitude modulation due to respiration and predict fluid responsiveness are disclosed. Disclosed embodiments include a method for assessing fluid responsiveness implemented in a digital computer with one or more processors comprising: (a) measuring a cardiovascular signal, and (b) computing a dynamic index predictive of fluid responsiveness from said cardiovascular signal using a nonlinear state space estimator. According to one particular embodiment, and without limitation, the nonlinear state space estimator is based on a model for cardiovascular signals such as arterial blood pressure or plethysmogram signals, and employs a marginalized particle filter to estimate a dynamic index predictive of fluid responsiveness that is substantially equivalent to a variation in pulse pressure of said cardiovascular signal.

CROSS-REFERENCE TO RELATED APPLICATIONS

This application claims the benefit of U.S. Provisional Application No.61/246,454 and U.S. Provisional Application No. 61/246,456 both filed on2009-09-28 by the present inventors, which are incorporated herein byreference.

TECHNICAL FIELD

Disclosed embodiments relate to methods and apparatuses for cardiacmonitoring. Specifically, they relate to methods and apparatuses for todynamic estimation of fluid responsiveness.

BACKGROUND

Indicators and methods for noninvasive determination of fluid status ofpatients are important for monitoring of the condition of critical carepatients. In many critical care settings clinicians must decide whetherpatients should be given intravenous fluid boluses and other therapiesto improve perfusion. Excessive fluid can be damaging by impairing lungfunction when it decreases oxygen delivery to tissues and contributes toorgan failure. Insufficient fluid can result in insufficient tissueperfusion which can also contribute to organ failure. Determining thebest course of fluid therapy for a given patient is difficult andclinicians have few clinical signs to guide them.

Fluid administration in hemodynamically unstable patients is often amajor challenge when it-comes to measuring stroke volume (SV), cardiacoutput (CO), or other hemodynamic parameters in real time. Correctclinical assessment of hypovolemia is difficult, as is the decision toundertake fluid resuscitation as the initial treatment strategy.Specifically, it is very difficult to predict whether a hemodynamicallyunstable patient will respond to fluid therapy with an increase instroke volume and cardiac output. Moreover, fluid overload can causesignificant pulmonary or cardiac dysfunction, whereas fluidinsufficiency may cause tissue damage resulting in vital organdysfunction. A patient's fluid responsiveness is the major and mostimportant determinant to assess the adequacy of fluid resuscitationtherapy and to ensure optimal cardiac performance and organ perfusion.

There are several dynamic parameters that can be used to assess fluidresponsiveness from arterial blood pressure (ABP) and in some cases fromplethysmogram signals. Several bedside indicators of ventricular preloadhave been used as predictors of fluid responsiveness. Right arterialpressure (RAP) and pulmonary artery occlusion pressure (PAOP) arecommonly used in the intensive care unit (ICU) when deciding toadminister fluids. Other bedside indicators of ventricular preloadinclude right ventricular end diastolic volume (RVEDV) and leftventricular end diastolic area (LVEDA) measured with transesophagealechocardiography. Several studies and case reports have shown, however,that these static indicators based on cardiac filling pressures manyhave poor predictive value and often fail to give adequate informationabout fluid responsiveness.

Several studies have confirmed the clinical significance of monitoringthe variations observed in left ventricular stroke volume that resultfrom the interaction of the cardiovascular system and the lungs undermechanical ventilation. These stroke volume variations (SVV) are causedby the cyclic increases and decreases in the intrathoracic pressure dueto the mechanical ventilation, which lead to variations in the cardiacpreload and afterload. SVV has recently been extensively investigatedand several studies have shown the usefulness of using SVV as predictorof fluid responsiveness in various clinical situations. Several otherparameters based on SVV have been found to be useful as well. Inparticular, systolic pressure variation (SPV) with its delta-Up anddelta-Down components has been found to be a very useful predictor offluid responsiveness. SPV is based on the changes in the arterial pulsepressure due to respiration-induced variations in stroke volume. Yetanother parameter that has recently been investigated and shown to be avalid indicator of fluid responsiveness is the pulse pressure variation(PPV).

The pulse pressure variation (PPV) index is a measure of the respiratoryeffect on the variation of systemic arterial blood pressure in patientsreceiving full mechanical ventilation. It is a dynamic predictor ofincreases in cardiac output due to an infusion of fluid. Numerousstudies have demonstrated that PPV is one of the most sensitive andspecific predictors of fluid responsiveness. Specifically, PPV has beenshown to be useful as a dynamic indicator to guide fluid therapy indifferent patient populations receiving mechanical ventilation. Forinstance, PPV was found to exhibit better performance as a predictor offluid responsiveness in patients before off-pump coronary artery bypassgrafting than standard static pre-load indices. PPV has also been shownto be useful for predicting and assessing the hemodynamic effects ofvolume expansion and a reliable predictor of fluid responsiveness inmechanically ventilated patients with acute circulatory failure relatedto sepsis.

The standard method for calculating PPV often requires simultaneousrecording of arterial and airway pressure. Pulse pressure (PP) ismanually calculated on a beat-to-beat basis as the difference betweensystolic and diastolic arterial pressure. Maximal PP (PP_(max)) andminimal PP (PP_(min)) are calculated over a single respiratory cycle,which is determined from the airway pressure signal. Pulse pressurevariations PPV are calculated in terms of PP_(max) and PP_(min), andexpressed as a percentage,

$\begin{matrix}{{P\; P\;{V(\%)}} = {100 \times \frac{{PP}_{\max} - {PP}_{\min}}{\left( {{PP}_{\max} + {PP}_{\min}} \right)/2}}} & (1)\end{matrix}$

There are few available methods or apparatuses to automatically estimatePPV accurately and reliably. A limitation of current PPV methods is thatthey do not work adequately in regions of abrupt hemodynamic changes andthis may limit their applicability in operating room environments.Additionally, they do not work for spontaneously breathing patients.

SUMMARY

Disclosed embodiments include a method for assessing fluidresponsiveness implemented in a digital computer with one or moreprocessors comprising: (a) measuring a cardiovascular signal, and (b)computing a dynamic index predictive of fluid responsiveness from saidcardiovascular signal using a nonlinear state space estimator. Accordingto one particular embodiment, and without limitation, the nonlinearstate space estimator is based on a model for cardiovascular signalssuch as arterial blood pressure or plethysmogram signals, and uses amarginalized particle filter to estimate a dynamic index predictive offluid responsiveness that is substantially equivalent to a variation inpulse pressure of said cardiovascular signal.

BRIEF DESCRIPTION OF THE DRAWINGS

Disclosed embodiments are illustrated by way of example, and not by wayof limitation, in the figures of the accompanying drawings.

FIG. 1 shows a block diagram of the method according to one embodiment.

FIG. 2 shows an estimation example based on the state space model andestimator according to one embodiment.

FIG. 3 shows an estimation example of PPV based on the state space modeland estimator according to one embodiment.

FIG. 4 shows a spectrogram of a sample signal used in the assessmentdataset.

FIG. 5 shows the performance results of the proposed method for PPVestimation according to one embodiment.

FIG. 6 illustrates the process of automatic PPV estimation duringartifacts according to one embodiment.

FIG. 7 illustrates the process of automatic PPV estimation duringartifacts according to one embodiment.

DETAILED DESCRIPTION

Disclosed embodiments include a method for assessing fluidresponsiveness 102 implemented in a digital computer with one or moreprocessors, said method comprising: (a) measuring a cardiovascularsignal 104, and (b) computing a dynamic index predictive of fluidresponsiveness from said cardiovascular signal using a nonlinear statespace estimator 106. The cardiovascular signal can be arterial bloodpressure (ABP), the photoplethysmographic signal (PLETH), central venouspressure (CVP), intracranial pressure (ICP), and others. According toone particular embodiment, and without limitation, the nonlinear statespace estimator is based on a model for cardiovascular signals 108 suchas arterial blood pressure or plethysmogram signals, and uses amarginalized particle filter 110 to estimate a dynamic index predictiveof fluid responsiveness 112 that is substantially equivalent to avariation in pulse pressure of said cardiovascular signal. The methodcan be implemented as part of an medical monitoring apparatus forassessing fluid responsiveness, comprising: (a) a measuring unit (i.e.physiological signal acquisition sensors and analog to digitalconverters) for measuring a cardiovascular signal 104; and (b) aprocessing unit for computing a dynamic index predictive of fluidresponsiveness from said cardiovascular signal using a nonlinear statespace estimator 106. According to a specific embodiment, the method canbe implemented in a medical system with one or more processors,physiological signal acquisition, analog to digital converters, memory,and output displays such as the typical bedside monitors used inclinical care. Alternatively, it can be implemented in a digitalcomputer with one or more processors to analyze physiological signalsand display or output, an indicator of fluid responsiveness that can beused to guide fluid therapy.

The following sections describe an embodiment of the a method forassessing fluid responsiveness 102 for the case of arterial bloodpressure signals (ABP) and estimation of the pulse pressure variationindex (PPV). Substantially equivalent embodiments without departing fromthe spirit of the disclosed embodiments are applicable to othercardiovascular signals and other derived parameters in addition to PPV.For instance, the state space model 108 and estimator 110 are applicableto the case of plethysmogram signals, and any other dynamic indexpredictive of fluid responsiveness substantially equivalent to avariation in pulse pressure of said cardiovascular signal due torespiration.

Certain specific details are set forth in the following description andfigures to provide a thorough understanding of various embodimentsdisclosed. Certain well-known details often associated with computingand software technology are not set forth in the following disclosure toavoid unnecessarily obscuring the various disclosed embodiments.Further, those of ordinary skill in the relevant art will understandthat they can practice other embodiments without one or more of thedetails described below. Aspects of the disclosed embodiments may beimplemented in the general context of computer-executable instructions,such as program modules, being executed by a computer, computer server,or device containing a processor. Generally, program modules includeroutines, programs, objects, components, data structures, etc. thatperform particular tasks or implement particular abstract data types.Aspects of the disclosed embodiments may also be practiced indistributed computing environments where tasks are performed by remoteprocessing devices that are linked through a communications network. Ina distributed computing environment, program modules may be located inboth local and remote storage media including memory storage devices.Those skilled in the art will appreciate that, given the description ofthe modules comprising the disclosed embodiments provided in thisspecification, it is a routine matter to provide working systems whichwill work on a variety of known and commonly available technologiescapable of incorporating the features described herein.

A. Mathematical Notation

In the detailed description of this application the followingmathematical notation is used. The mathematical notation is thoroughlydescribed in this section to aid those skilled in the art to understand,implement, and practice the disclosed embodiments. We use boldface todenote random processes, normal face for deterministic parameters andfunctions, upper case letters for matrices, lower case letters forvectors and scalars, superscripts in parenthesis for particle indices,upper-case superscripts for nonlinear/linear indication, and subscriptsfor time indices. For example, the nonlinear portion of the state vectorfor the i^(th) state trajectory (i.e., particle) is denoted as χ_(n)^(N,(i)) where n represents the discrete time index and (i) denotes thei^(th) particle. The unnormalized particle weights are denoted as {tildeover (w)}^((i)) and the normalized particle weights as w^((i)). Thestate trajectories before resampling are denoted as {tilde over (χ)}_(n)^((i)) and as χ_(n) ^((i)) after resampling.

B. State-Space Model

According to one embodiment of the proposed method and apparatus forassessing fluid responsiveness 102, and without limitation, themethod/apparatus utilizes a sequential Monte Carlo estimation method 110which is based on the state-space model approach. The state-space modelis a mathematical expression to describe the evolution of any physicalsystem's unobservable state χ_(n) and its relation to measurement y_(n),where the state χ_(n) is a vector of parameters representing thesystem's condition. The state-space method is a technique to estimatethe state χ_(n) as a function of measurement y_(n) utilizing thestate-space model. The typical state-space model can be expressed as,χ_(n+1)ƒ(χ_(n))+u _(n)  (2)y _(n) =h(χ_(n))+v _(n)  (3)where (2) is a process model, (3) a measurement model, ƒ (•) and h (•)functions of the state χ_(n), and u_(n) and v_(n) uncorrelated whitenoise with variances q and r. A designer needs to incorporate priordomain knowledge of a system into the state-space model and define thefunctions ƒ (•) and h (•). The flexibility and versatility of thestate-space method are attributable to two functions, which can beeither linear or nonlinear.B.1. Measurement Model

According to one embodiment of the proposed method and apparatus forassessing fluid responsiveness 102, and without limitation, themeasurement model of the cardiovascular signal is shown in (4)-(6),where γ_(n) is the respiratory signal, μ_(n), the amplitude-modulatedcardiac signal, ρ_(k,n) the amplitude modulation factor of the k^(th)cardiac harmonic partial, κ_(k,n) the k^(th) cardiac harmonic partial,θ_(n) ^(r) the instantaneous respiratory angle, θ_(n) ^(c) theinstantaneous cardiac angle, ƒ_(n) ^(r) the instantaneous respiratoryrate, N_(h) ^(r) the number of respiratory partials, N_(h) ^(c) thenumber of cardiac partials, and v_(n) the white Gaussian measurementnoise with variance r.

$\begin{matrix}{y_{n} = {{\gamma_{n} + \mu_{n} + v_{n}} = {\gamma_{n} + {\sum\limits_{k = 1}^{N_{h}^{c}}{\rho_{k,n}\kappa_{k,n}}} + v_{n}}}} & (4) \\{\gamma_{n} = {{\sum\limits_{k = 1}^{N_{h}^{r}}{r_{1,k,n}{\cos\left( {k\;\theta_{n}^{r}} \right)}}} + {r_{2,k,n}{\sin\left( {k\;\theta_{n}^{r}} \right)}}}} & (5) \\{\rho_{k,n} = {1 + {\sum\limits_{j = 1}^{N_{h}^{r}}{m_{1,k,j,n}{\cos\left( {j\;\theta_{n}^{r}} \right)}}} + {m_{2,k,j,n}{\sin\left( {j\;\theta_{n}^{r}} \right)}}}} & (6) \\{\kappa_{k,n} = {{\sum\limits_{k = 1}^{N_{h}^{c}}{c_{1,k,n}{\cos\left( {k\;\theta_{n}^{c}} \right)}}} + {c_{2,k,n}{\sin\left( {k\;\theta_{n}^{c}} \right)}}}} & (7)\end{matrix}$B.2. Process Model

According to one embodiment of the proposed method and apparatus forassessing fluid responsiveness 102, and without limitation, the processmodel describes the evolution of each element of the state χ_(n). In ourapplication, χ_(n) includes the instantaneous respiratory and cardiacangles θ_(n) ^(r) and θ_(n) ^(c), the instantaneous mean cardiacfrequency ƒ _(n) ^(c), the instantaneous respiratory and cardiacfrequencies ƒ_(n) ^(r) and ƒ_(n) ^(c), and the sinusoidal coefficients{r_(1,k,n), r_(2,k,n), c_(1,k,n), c_(2,k,n), m_(1,k,n), m_(2,k,n)}. Theprocess model can be expressed as,θ_(n+1) ^(r)=2π(n+1)T _(s)ƒ^(r)  (8)θ_(n+1) ^(c)=θ_(n) ^(c)+2πT _(s)ƒ_(n) ^(c)  (9)ƒ _(n+1) ^(c) =g[ ƒ _(n) ^(c) +u _(ƒ) _(c) _(,n)]  (10)ƒ_(n+1) ^(c)= ƒ _(n) ^(c)+α(ƒ_(n) ^(c)− ƒ _(n) ^(c))+u _(ƒ) _(c)_(,n)  (11)r. _(,k,n+1) =r. _(,k,n) +u _(r,n)  (12)c. _(,k,n+1) =c. _(,k,n) +u _(c,n)  (13)m. _(,k,n+1) =m. _(,k,n) +u _(m,n)  (14)where ƒ^(r) is the constant respiratory frequency, ƒ_(n) ^(c) theinstantaneous cardiac frequency, T_(s) the sampling period, ƒ _(n) ^(c)the instantaneous mean cardiac frequency, α the autoregressivecoefficient, and u_(r,n), u_(c,n), and u_(m,n), the process noises withvariances q_(r), q_(c), and q_(m). Since this process model is for ABPsignals of patients with full mechanical support, the respiratoryfrequency ƒ^(r) has a known constant value. The clipping function g[•]limits the range of instantaneous mean frequencies, which can be writtenas,

$\begin{matrix}{{g\lbrack f\rbrack} = \left\{ \begin{matrix}{f_{\max} - \left( {f - f_{\max}} \right)} & {f_{\max} < f} \\f & {f_{\min} < f \leq f_{\max}} \\{f_{\min} + \left( {f_{\min} - f} \right)} & {f \leq {f_{\min}.}}\end{matrix} \right.} & (15)\end{matrix}$

The range of instantaneous mean frequencies is assumed to be known asdomain knowledge. Other functions could be used here that are generallysigmoidal in shape and limit the range of instantaneous frequencies suchas arctangent and hyperbolic tangent functions.

C. State Estimation

According to one embodiment of the proposed method and apparatus forassessing fluid responsiveness 102, and without limitation, thenonlinear state space estimator is a sequential Monte Carlo method(SMCM). The SMCM is a computational Bayesian approach to estimation ofan unknown state χ_(n), recursively given measurement y_(n) and thestate-space model. The SMCM estimates the posterior mean of the state E[p(χ_(0:n)|y_(0:n))] as a weighted sum of random samples χ_(0:n) ^((i)),which can be expressed as,

$\begin{matrix}{{E\left\lbrack {p\left( {x_{0:n}❘y_{0:n}} \right)} \right\rbrack} = {\sum\limits_{i = 1}^{N_{p}}{x_{0:n}^{(i)}w_{0:n}^{(i)}}}} & (16)\end{matrix}$where χ_(0:n)={χ₀, . . . , χ_(n),} and w_(0:n) ^((i)) is the randomsample weight. The SMCM is better known as a particle filter (PF), sincethe random samples are often referred to as particles. The proposed SMCMovercomes the limitations of the conventional SMCM such as samplingdegeneracy and impoverishment issues. It also addresses properly themulti-modality issue of the posterior distribution p(ƒ_(n)|y_(0:n))where ƒ_(n) is the fundamental frequency of a multiharmonic signa. Theproposed SMCM is called, maximum a-posteriori adaptive marginalizedparticle filter (MAM-PF). The method is described in detail in Table(Algorithm 1).

The MAM-PF ABP tracker continuously estimates the parameters of theprocess model shown in (8)-(13). The number of parameters to track is4N_(h) ^(r)+2N_(c) ^(c)+1. Obtaining accurate estimation of thoseparameters is the first step toward PPV index estimation.

D. ABP and Pleth Signal Envelope Estimation

Given the estimated signal parameters in (8)-(14), it is possible toestimate the upper envelope (e_(μ,n)) and lower envelope (e_(l,n)) ofABP or pleth signals by following steps below,

$\theta_{\max}^{c} = {\arg{\max\limits_{\theta}{\sum\limits_{k = 1}^{N_{h}^{c}}{\rho_{k,n}\left\lbrack {{c_{1,k,n}{\cos\left( {k\;\theta} \right)}} + {c_{2,k,n}{\sin\left( {k\;\theta} \right)}}} \right\rbrack}}}}$$\theta_{\min,n}^{c} = {\arg{\min\limits_{\theta}{\sum\limits_{k = 1}^{N_{h}^{c}}{\rho_{k,n}\left\lbrack {{c_{1,k,n}{\cos\left( {k\;\theta} \right)}} + {c_{2,k,n}{\sin\left( {k\;\theta} \right)}}} \right\rbrack}}}}$κ_(max , k, n) = c_(1, k, n)cos (k θ_(max , n)^(c)) + c_(2, k, n)sin (k θ_(max , n)^(c))κ_(min , k, n) = c_(1, k, n)cos (k θ_(min , n)^(c)) + c_(2, k, n)sin (k θ_(min , n)^(c))$\begin{matrix}{e_{\mu,n} = {\gamma_{n} + {\sum\limits_{k = 1}^{N_{h}^{c}}{\rho_{k,n}\kappa_{\max,k,n}}}}} & (17) \\{e_{l,n} = {\gamma_{n} + {\sum\limits_{k = 1}^{N_{h}^{c}}{\rho_{k,n}\kappa_{\min,k,n}}}}} & (18)\end{matrix}$where arg max_(χ)ƒ(χ) and arg min_(χ)ƒ(χ) are operators to obtain thevalue of χ for which ƒ(χ) attains its maximum and minimum values,respectively. The top plot in FIG. 2 shows a five respiratory cycleperiod of an ABP signal y_(n) its estimate ŷ_(n), and its estimatedenvelopes e_(μ,n) and e_(l,n).

Algorithm 1 Table-MAM-PF Method Steps Initialization  for  i = 1, …  , N_(p)  do     Sample  x₀^(N, (i)) ∼ π₀(x₀^(N))    x̂_(0 : −1)^(L, (i)) = E[x₀^(L, (i))|x₀^(N, (i))]   end  for  for  i = 1, …  , N_(p)  do    α₀^((i)) = π₀(x₀^(N, (i)))p(y₀|x₀^(N, (i)), x₀^(L, (i)))    z₀^((i)) = x₀^((i))   end  for$\mspace{25mu}{i^{*} = {{{{\underset{i}{argmax}\mspace{14mu}\alpha_{0}^{(i)}}\&}\mspace{14mu}{\hat{x}}_{0}} = x_{0}^{i^{*}}}}$  for  n = 1, …  , N_(T)  do     for  i = 1, …  , N_(p)  do     MAP  Estimation      Particle  Propagation       x_(n)^(N, (i)) ∼ q_(n)(x_(n)^(N, (i))|x_(n − 1)^(N, (i)), y_(n))$\mspace{101mu}{k^{*} = {\underset{k}{argmax}\mspace{14mu}\alpha_{n - 1}^{(k)}{p\left( x_{n}^{N,{(i)}} \middle| x_{n - 1}^{N,{(k)}} \right)}}}$     Marginalized  Sequential  Estimation      Measurement  Update       ŷ_(n|0 : n − 1) = H_(n)(x_(n)^(N, (i)))x̂_(n|0 : n − 1)^(L, (k^(*)))       e_(n) = y_(n) − ŷ_(n|0 : n − 1)       R_(v, n) = [e_(n)e_(n)^(T) − H_(n)(x_(n)^(N, (i)))C_(n|0 : n − 1)^((k^(*)))H_(n)(x_(n)^(N, (i)))^(T)]₊       R̂_(v, n) = β R̂_(v, n − 1)^((k^(*))) + (1 − β)R_(v, n)       R_(e, n) = H_(n)(x_(n)^(N, (i)))C_(n|0 : n − 1)^((k^(*)))H_(n)(x_(n)^(N, (i)))^(T) + R̂_(v, n)       K_(n) = C_(n|0 : n − 1)^((k^(*)))H_(n)(x_(n)^(N, (i)))^(T)(R_(e, n))⁻¹       x̂_(n|0 : n)^(L) = x̂_(n|0 : n − 1)^(L, (k^(*))) + K_(n)e_(n)       C_(n|0 : n) = [I − K_(n)H_(n)(x_(n)^(N, (i)))]C_(n|0 : n − 1)^((k^(*)))     Time  Update        R̂_(v, n)^((i)) = R̂_(v, n)       x̂_(n + 1|0 : n)^(L, (i)) = F_(n)(x_(n)^(N, (i)))x̂_(n|0 : n)^(L)       C_(n + 1|0 : n)^((i)) = F_(n)(x_(n)^(N, (i)))C_(n|0 : n)F_(n)(x_(n)^(N, (i)))^(T) + Q_(u)^(L)     α_(n)^((i)) = α_(n − 1)^((k^(*)))p(y_(n)|x_(n)^(N, (i)), x̂_(n|0 : n − 1)^(L, (k^(*))))p(x_(n)^(N, (i))|x_(n − 1)^(N, (k^(*))))     x̂_(n)^((i)) = [x̂_(n|0 : n)^(L), x_(n)^(N, (i))]^(T)     z_(0 : n)^((i)) = [z_(0 : n − 1)^((k^(*))), x̂_(n)^((i))]   end  for  Update  MAP  State  Estimate$\mspace{25mu}{i^{*} = {\underset{i}{argmax}\alpha_{n}^{(i)}}}$  x̂_(0 : n) = z_(0 : n)^(i^(*)) end  forE. Pulse Pressure Signal Envelope Estimation

The pulse pressure (PP) signal is the difference between the upperenvelope e_(μ,n) and lower envelope e_(l,n) of the cardiovascularpressure signal (ABP, ICP, pleth, ECG, etc). This PP signal oscillatesroughly at the respiratory rate as shown in the bottom plot in FIG. 2.This oscillation is due to the respiratory effect on the variation ofsystemic ABP under full ventilation support. Within each respiratorycycle PP reaches its maximum (PP_(max)) and minimum (PP_(min)) values,which are two critical parameters to compute the PPV index.

Traditionally, the PP_(max) and PP_(min), values have been computed onlyonce per each respiratory cycle. Given the estimated signal parametersin (8)-(14), however, according to one embodiment, the continuousequivalents of PP_(max) and PP_(min) are computed. They are the upperε_(μ,n) and lower ε_(l,n) envelopes of the PP signal. The upper envelopeis the continuous estimate of PP_(max) and the lower envelope ε_(l,n)that of PP_(min). The ε_(μ,n) and ε_(l,n) values can be estimated asdescribed below,

$\varrho_{k,n} = {{\sum\limits_{j = 1}^{N_{h}^{r}}{m_{1,k,j,n}{\cos\left( {j\;\theta} \right)}}} + {m_{2,k,j,n}{\sin\left( {j\;\theta} \right)}}}$$\theta_{\max,n}^{r} = {\arg\;{\max\limits_{\theta}{\sum\limits_{k = 1}^{N_{h}^{c}}{\left( {1 + \varrho_{k,n}} \right)\left( {\kappa_{\max,k,n} - \kappa_{\min,k,n}} \right)}}}}$$\theta_{\min,n}^{r} = {\arg\;{\min\limits_{\theta}{\sum\limits_{k = 1}^{N_{h}^{c}}{\left( {1 + \varrho_{k,n}} \right)\left( {\kappa_{\max,k,n} - \kappa_{\min,k,n}} \right)}}}}$$\varrho_{\max,k,n} = {{\sum\limits_{j = 1}^{N_{h}^{r}}{m_{1,k,j,n}{\cos\left( {j\;\theta_{\max,n}^{r}} \right)}}} + {m_{2,k,j,n}{\sin\left( {j\;\theta_{\max,n}^{r}} \right)}}}$$\varrho_{\min,k,n} = {{\sum\limits_{j = 1}^{N_{h}^{r}}{m_{1,k,j,n}{\cos\left( {j\;\theta_{\min,n}^{r}} \right)}}} + {m_{2,k,j,n}{\sin\left( {j\;\theta_{\min,n}^{r}} \right)}}}$$\begin{matrix}{ɛ_{\mu,n} = {\sum\limits_{k = 1}^{N_{h}^{c}}{\left( {1 + \varrho_{\max,k,n}} \right)\left( {\kappa_{\max,k,n} - \kappa_{\min,k,n}} \right)}}} & (19) \\{ɛ_{l,n} = {\sum\limits_{k = 1}^{N_{h}^{c}}{\left( {1 + \varrho_{\min,k,n}} \right)\left( {\kappa_{\max,k,n} - \kappa_{\min,k,n}} \right)}}} & (20)\end{matrix}$where 1+

is equal to ρ_(k,n), and ε_(μ,n) and ε_(l,n) are the continuousestimates of the PP_(max) and PP_(min), respectively. The lines in thebottom plot in FIG. 2 represent the upper ε_(μ,n) and lower ε_(l,n)envelopes of the PP signal obtained by following the method describedabove.F. Pulse Pressure Variation Calculation

According to one embodiment of the proposed method and apparatus forassessing fluid responsiveness 102, and without limitation, given theε_(μ,n) and lower ε_(l,n) values, the PPV index is calculated directlyas follows,

$\begin{matrix}{{P\; P\;{V(\%)}} = {100 \times \frac{ɛ_{\max} - ɛ_{\min}}{\left( {ɛ_{\max} + {ɛ_{\min}/2}} \right)}}} & (21)\end{matrix}$

According to one embodiment, PPV is typically calculated over severalrespiratory cycles (at least one cycle is required). This new PPV indexis different from the traditional PPV index described in (1) because thenew one is continuous in time while the traditional one can be obtainedonly once per each respiratory cycle.

FIG. 3 illustrates an example of the automatically computed continuousPPV index and the manually obtained discrete PPV index of a real 10 minABP signal. Each hollow white dot represents a “discrete” PPV index,which can be obtained once per each respiratory cycle.

G. Validation and Testing of the Method According to One Embodiment.

The following description presents the experimental results and tests ofa prototype according to a disclosed embodiment.

G.1. Method Validation: Assessment Data The Massachusetts GeneralHospital (MGH) waveform database on PhysioNet is a comprehensivecollection of electronic recordings of hemodynamic andelectrocardiographic waveforms patients in critical care units. Itconsists of recordings from 250 patients representing

TABLE 1 Summary of user-specified design parameters for the PPV indexestimator Name Symbol Value No. particles N^(p) 250  No. cardiaccomponents N^(c) 10 No. respiratory components N^(r)  3 Minimum heartrate f_(min) ^(c)  50/60 Hz Maximum heart rate f_(max) ^(c) 150/60 HzMeasurement noise variance r var(y)/1e3 Cardiac frequency varianceq_(f ) _(c) 7e−6 T_(s) Respiratory amplitude variance q_(a)&q_(b) 1e−5T_(s) Modulation factor amplitude variance q_(c)&q_(d) 1e−9 T_(s)Cardiac amplitude variance q_(e)&q_(f) 1e−6 T_(s) Initial respiratoryamplitude u_(a)&u_(b) std(y)/1e2 Initial modulation factor amplitudeu_(c)&u_(d) std(y)/1e2 Initial cardiac amplitude u_(e)&u_(f) std(y)a broad spectrum of physiologic and pathophysiologic states. The typicalrecording includes arterial blood pressure (ABP) signals in addition toseven other types of waveforms. By visually inspecting the spectrogramof ABP signals we identified 26 patients whose respiratory rate remainedconstant at least for 10 consecutive minutes. We used the constantrespiratory rate shown in spectrogram as an indicator of fullrespiratory support. FIG. 4 shows the spectrogram of one of the 26 ABPsignals. The total duration of the 26 ABP signals was a little bit over4 hours (260 min). The original sample rate ƒ_(s) of the signals was 360Hz, but they were downsampled by a factor of 9, so that the final samplerate ƒ_(s) was 40 Hz.

Typically the constant respiratory rate ƒ^(r) of subjects under fullventilation support is easily accessible medical information since it isequal to the machine ventilation rate.

In this embodiment, the number of cardiac components N^(c) was 10 andthat of respiratory components N^(r) 3. The number of particles was 250.Table 1 lists the parameter values used for the PPV index estimator.

G.2. Method Validation: Manual PPV Annotations (Gold Standard)

We manually detected the peaks and troughs of the ABP signals andcalculated the PPV indices (gold standard) as defined in (1). They arereferred to as manual PPV indices PPV_(manu). PPV_(auto) represents PPVindices obtained using the proposed PPV index tracking algorithm.

We had to exclude 3 ABP signals out of the 26 ABP signals from ourvalidation study because it was not possible to annotate either those 3ABP signals or their PP signals due to a high noise level or abnormalcardiac activity.

G. 3. Method Validation: Statistical Analysis

We took 5 PPV index measurements per subject separately by 2 min. EachPPV index measurement is an averaged value over 5 respiratory cycles.FIG. 3 shows the 2 min apart measurement periods. The proposed PPV indextracking algorithm was assessed by calculating the agreement(mean±standard deviation) between PPV_(auto) and PPV_(manu) measurementsand using Bland-Altman analysis.

A Bland-Altman plot is a statistical and visualization method that isoften used in the assessment of PPV estimation method in order todetermine the agreement between two different PPV estimates across thePPV range. It shows the difference ΔPPV between PPV_(auto) andPPV_(manu) on the y-axis and the PPV_(manu) on the x-axis. It visualizesnot only the overall accuracy of estimation but also any estimation biasor trend versus PPV_(manu). We used it to compare the gold standardusing manual annotations with our automatic estimation algorithm.

G. 4. Method Validation: Results

FIG. 5 depicts the Bland-Altman plot of the 23 subjects. There are 5 PPVmeasurements available per each subject. In this embodiment, most (98%)of PPV_(auto) measurements were in agreement with PPV_(manu)measurements within ±5% accuracy, and 100% within ±8% accuracy.

TABLE 2 Summary of the mean and standard deviation of the PPV_(manu) andPPV_(auto) Subject PPV_(manu) (%) PPV_(auto) (%) 1 23.2 ± 1.0 23.6 ± 1.02 23.0 ± 0.3 22.8 ± 1.6 3 13.7 ± 0.8 14.3 ± 1.2 4 12.3 ± 1.1  9.1 ± 3.15 36.3 ± 5.9 33.9 ± 6.7 6  8.7 ± 1.8  9.5 ± 2.0 7 18.2 ± 2.6 18.8 ± 1.78  6.3 ± 1.3  6.2 ± 0.8 9  9.6 ± 1.7  7.1 ± 0.8 10 15.9 ± 0.8 14.8 ± 0.811  5.3 ± 1.0  5.5 ± 0.7 12  3.1 ± 0.2  4.0 ± 0.3 13 31.8 ± 1.6 31.5 ±0.8 14 15.3 ± 2.8 13.1 ± 1.8 15 33.9 ± 1.1 34.3 ± 1.8 16  4.6 ± 1.1  3.4± 1.0 17 20.4 ± 2.3 20.1 ± 2.7 18  9.8 ± 0.8 11.9 ± 2.2 19  9.6 ± 2.111.2 ± 1.4 20 11.9 ± 0.5 10.3 ± 1.2 21  8.7 ± 1.8 11.9 ± 1.9 22 21.7 ±1.9 18.9 ± 1.2 23 20.9 ± 2.1 20.0 ± 0.7Table 2 lists the mean±standard deviation of 5 PPV_(manu) and PPV_(auto)measurements for each subject. The second column is for PPV_(manu) andthe third column for PPV_(auto).

A design difference of the proposed method with respect to previouslypublished method is the fact that the proposed method is based on astatistical state-space model and estimation of the cardiovascularpressure pressure signal. The state-space modeling stage results in analgorithm that is more robust to hemodynamic changes and artifacts. Thestatistical state-space signal model and associated model parameterestimation algorithm automatically filter out noise and artifact thatcannot be captured with the model. Since the statistical signal model isbased on cardiovascular physiology and pathophysiology, signal featuresthat are not physiological in nature are automatically filtered out.Additionally, the model is general enough to accurately model botharterial blood pressure signals and plethysmogram signals. Consequently,it can also be used to calculate the pleth variability index (PVI).

FIG. 6 and FIG. 7 exemplify a case where signal features that are notphysiological in nature are automatically filtered out resulting in moreaccurate PPV index estimation than manual annotation. The top plot inFIG. 6 illustrates 6 respiratory cycles of the ABP signal of Subject 3and its estimate. It also shows the manually annotated signal envelopesand the automated computed signal envelopes. The bottom plot in FIG. 7depicts the PPV_(manu) and PPV_(auto) over the same period. Around 399s, the PPV_(manu) value abruptly increases up to 78% while thePPV_(auto) value remains at 20%. Around 403 s, the PPV_(manu) valuereturns to 20%. FIG. 7 focuses on the time period marked with the blackrectangular box in FIG. 6. The top plot in FIG. 7 shows that the heartbeat between 398 s and 399 s is irregular and abnormal. As a result, thecorresponding PP_(manu) shown in the bottom plot reaches a very smallminimum value (PP_(min,manu): 38%) between 398 s and 399 s. However, theautomatically computed minimum PP value (PP_(min,auto)) over the sameperiod is as high as 70%. This discrepancy between the manual annotationand the proposed automatic method results from the capability of theMAM-PF algorithm, which estimates the ABP signal based on thestate-space model. While the original heart beat 398 s and 399 s in FIG.6 is abnormal in a physiological sense, the estimated heart beat overthe same time period shows the physiologically expected morphology andlocation of the heart beat.

H. Method for Assessing Fluid Responsiveness in Spontaneously BreathingPatients

According to one embodiment, the method and apparatus for assessingfluid responsiveness is adapted to work for spontaneously breathingpatients (i.e. patients not mechanically ventilated). In this section wedescribe an embodiment for spontaneously breathing patients. Fluidresponsiveness prediction in spontaneously breathing subjects isdifficult. Since the swings in intrathoracic pressure are minor duringspontaneous breathing, dynamic parameters like pulse PPV and systolicpressure variation (SPV) are usually small and more difficult tocalculate.

H.1. Measurement Model

According to one embodiment, the measurement model of the ABP signal isshown in (22)(25), where γ_(n) is the respiratory signal, μ_(n) theamplitude-modulated cardiac signal, ρ_(k,n) the amplitude modulationfactor of the k^(th) cardiac harmonic partial, κ_(k,n) the k^(th)cardiac harmonic partial, θ_(n) ^(r) the instantaneous respiratoryangle, θ_(n) ^(c) the instantaneous cardiac angle, ƒ_(n) ^(r) theinstantaneous respiratory rate, N_(h) ^(r) the number of respiratorypartials, N_(h) ^(c) the number of cardiac partials, and v_(n) the whiteGaussian measurement noise with variance r.

$\begin{matrix}{y_{n} = {{\gamma_{n} + \mu_{n} + \upsilon_{n}} = {\gamma_{n} + {\sum\limits_{k = 1}^{N_{h}^{c}}{\rho_{k,n}\kappa_{k,n}}} + \upsilon_{n}}}} & (22) \\{\gamma_{n} = {{\sum\limits_{k = 1}^{N_{h}^{r}}{r_{1,k,n}{\cos\left( {k\;\theta_{n}^{r}} \right)}}} + {r_{2,k,n}{\sin\left( {k\;\theta_{n}^{r}} \right)}}}} & (23) \\{\rho_{k,n} = {1 + {\sum\limits_{j = 1}^{N_{h}^{r}}{m_{1,k,j,n}{\cos\left( {j\;\theta_{n}^{r}} \right)}}} + {m_{2,k,j,n}{\sin\left( {j\;\theta_{n}^{r}} \right)}}}} & (24) \\{\kappa_{k,n} = {{\sum\limits_{k = 1}^{N_{h}^{c}}{c_{1,k,n}{\cos\left( {k\;\theta_{n}^{c}} \right)}}} + {c_{2,k,n}{\sin\left( {k\;\theta_{n}^{c}} \right)}}}} & (25)\end{matrix}$H.2. Process Model

The process model describes the evolution of each element of the stateχ_(n). In our application, χ_(n) includes the instantaneous respiratoryand cardiac angles θ_(n) ^(r) and θ_(n) ^(c), the instantaneous meanrespiratory ƒ _(n) ^(r) and cardiac ƒ _(n) ^(c) frequency, theinstantaneous respiratory ƒ_(n) ^(r) and cardiac ƒ_(n) ^(c) frequencies,and the sinusoidal coefficients {r_(1,k,n), r_(2,k,n), c_(1,k,n),c_(2,k,n), m_(1,k,n), m_(2,k,n)}. The process model can be expressed as,θ_(n+1) ^(r)=θ_(n) ^(r)+2πT _(s)ƒ_(n) ^(r)  (26)θ_(n+1) ^(c)=θ_(n) ^(c)+2πT _(s)ƒ_(n) ^(c)  (27)ƒ _(n+1) ^(r) =g[ ƒ _(n) ^(r) +u _(ƒ) _(r) _(,n)]  (28)ƒ _(n+1) ^(c) =g[ ƒ _(n) ^(c) +u _(ƒ) _(c) _(,n)]  (29)ƒ_(n+1) ^(r)= ƒ _(n) ^(r)+α(ƒ_(n) ^(r)− ƒ _(n) ^(r))+u _(ƒ) _(r)_(,n)  (30)ƒ_(n+1) ^(c)= ƒ _(n) ^(c)+α(ƒ_(n) ^(c)− ƒ _(n) ^(c))+u _(ƒ) _(c)_(,n)  (31)r. _(,k,n+1) =r. _(,k,n) +u _(r,n)  (32)c. _(,k,n+1) =c. _(,k,n) +u _(c,n)  (33)m. _(,k,n+1) =m. _(,k,n) +u _(m,n)  (34)where ƒ_(n) ^(r) is the instantaneous respiratory frequency, ƒ_(n) ^(c)the instantaneous cardiac frequency, T_(s) the sampling period, ƒ _(n)^(r) the instantaneous mean respiratory frequency, ƒ _(n) ^(c) theinstantaneous mean cardiac frequency, α the autoregressive coefficient,and u_(r,n), u_(c,n), and u_(m,n) the process noises with variancesq_(r), q_(c), and q_(m). The clipping function g[•] limits the range ofinstantaneous mean frequencies, which can be written as,

$\begin{matrix}{{g\lbrack f\rbrack} = \left\{ \begin{matrix}{f_{\max} - \left( {f - f_{\max}} \right)} & {f_{\max} < f} \\f & {f_{\min} < f \leq f_{\max}} \\{f_{\min} + \left( {f_{\min} - f} \right)} & {f \leq {f_{\min}.}}\end{matrix} \right.} & (35)\end{matrix}$The range of instantaneous mean frequencies is assumed to be known asdomain knowledge.H.3. Maximum A-Posteriori MPF

The proposed automated PPV index estimation method requires accurateestimates of the instantaneous respiratory frequency ƒ_(n) ^(r), theinstantaneous cardiac frequency ƒ_(n) ^(c), and the morphology of an ABPsignal. Within the state-space method framework, the Optimal MAM-PFcomputes the “optimal” trajectory of the state χ_(n). The Fast MAM-PF isan approximation of the Optimal MAM-PF, which requires dramatically lesscomputational burden. However, the Fast MAM-PF performs as well as theOptimal MAM-PF. Under full mechanical support, the respiratory rate ofsubjects is equal to the mechanical ventilation rate, which is known andconstant. Therefore, the Fast MAM-PF has to track only the instantaneouscardiac frequency ƒ_(n) ^(c) along with the signal morphology.

The ABP signal tracker has to track both the instantaneous respiratoryfrequency ƒ_(n) ^(r) and the instantaneous cardiac frequency ƒ_(n) ^(c)along with the signal morphology. Although the Fast MAM-PF based ABPsignal tracker is capable of tracking multiple frequencies, there aretwo major issues in using the Fast MAM-PF algorithm as the ABP signaltracker for ABP signals of spontaneously breathing subjects. The firstissue is that the morphology of the signal, which is represented by thesinusoidal coefficients in (24)-(25), does not belong to the linearstate any more. Since the modulating signal ρ_(k,n), is multiplied tothe cardiac signal κ_(k,n), their sinusoidal coefficients c._(,k,n) andm._(,k,j,n) are nonlinear parameters of the measurement model in (22).The Fast MAM-PF is applicable only to state-space models whose statevector can be partitioned into the linear and nonlinear portions. Thesecond issue is that the computational burden of the Fast MAM-PFincreases exponentially as the dimension of the state where particlefilters are used increases. The portion of the state space whereparticle filters are used is called the particle space. Since the newABP signal tracker has to estimate both the instantaneous respiratoryfrequency ƒ_(n) ^(r) and the instantaneous cardiac frequency ƒ_(n) ^(c),the dimension of the particle state becomes 2, which results in aquadruple increase of computational burden if the Fast MAM-PF has to beused for the current application. Nevertheless, it is a possibleembodiment that may be useful in certain situations. In order to addressthese two major issues we propose a new ABP signal tracker, which is amodified version of the Fast MAM-PF. It is called, the Dual MAM-PF. Theterm “Dual” reflects the fact that the particle space istwo-dimensional. While the Fast MAM-PF treats a two-dimensional particlespace as a whole, the Dual MAM-PF partitions the two-dimensionalparticle space into two one-dimensional particle spaces assumingindependence between two particle space variables, which are theinstantaneous respiratory frequency ƒ_(n) ^(r) and the instantaneouscardiac frequency ƒ_(n) ^(c.)

Suppose that the state vector χ can be partitioned as follows,

$\begin{matrix}{x_{n} = \begin{bmatrix}x_{n}^{P} \\x_{n}^{K}\end{bmatrix}} & (36)\end{matrix}$where χ_(n) ^(P) represents the particle state and χ_(n) ^(K) the Kalmanstate. The particle state is the portion of the state where particlefilters are used as defined earlier while the Kalman state is theportion of the state where extended Kalman filters are used. The statevariables whose posterior distributions are known to be multi-modalbelong to the particle state while those whose posterior distributionsare known to be Gaussian or uni-modal belong to the Kalman state. Giventhe state-space model in (4)-(7), instantaneous respiratory frequencyƒ_(n) ^(r) and the instantaneous cardiac frequency ƒ_(n) ^(c) are theparticle state variables and the sinusoidal coefficients such asr._(,k,n), c._(,k,n), and m._(,k,j,n) are the Kalman state variables.Assuming that the particle state variables are independent of each otherthe particle state χ_(n) ^(P) can be partitioned further as,

$\begin{matrix}{x_{n}^{P} = \begin{bmatrix}x_{n}^{P_{1}} \\x_{n}^{P_{2}}\end{bmatrix}} & (37) \\{x_{n + 1}^{P_{1}} = {f_{1,n}\left( {x_{n}^{P_{1}},u_{n}^{P_{1}}} \right)}} & (38) \\{x_{n + 1}^{P_{2}} = {f_{2,n}\left( {x_{n}^{P_{2}},u_{n}^{P_{2}}} \right)}} & (39)\end{matrix}$where χ_(n) ^(P) ¹ and χ_(n) ^(P) ² represent the first and secondparticle state variables, respectively. This partitioning breaks down atwo-dimensional particle space χ_(n) ^(P) into two one-dimensionalparticle spaces. The total posterior distribution is given by,

$\begin{matrix}\begin{matrix}{{p\left( {x_{0:n}❘y_{0:n}} \right)} = {{p\left( {{x_{0:n}^{K}❘y_{0:n}},x_{0:n}^{P}} \right)}{p\left( {x_{0:n}^{P}❘y_{0:n}} \right)}}} \\{= {{p\left( {{x_{0:n}^{K}❘y_{0:n}},x_{0:n}^{P}} \right)}{p\left( {x_{0:n}^{P_{1}}❘y_{0:n}} \right)}{{p\left( {x_{0:n}^{P_{2}}❘y_{0:n}} \right)}.}}}\end{matrix} & \begin{matrix}(40) \\(41)\end{matrix}\end{matrix}$

Algorithm 2 provides a complete description of the Dual MAM-PFalgorithm, where N_(T) represents the total number of signal samples,N_(p) the number of particles for each one-dimensional particle space, jthe particle state variable index, and i_(j) the particle index of thej^(th) particle state variable. The total number of particle used in theDual MAM-PF algorithm is 2N_(p) instead of N_(p) ². At each time step nthe Dual MAM-PF searches for the best trajectory of each particle i_(j)from the previous trajectory. This searching step can be written as,

$\begin{matrix}\begin{matrix}{k_{j}^{*} = {\underset{k_{j}}{argmax}\;\alpha_{j,{n - 1}}^{({k,_{j}})}{p\left( {{y_{n}❘x_{n}^{P_{j},{(i_{j})}}},{\hat{x}}_{n❘{0:{n - 1}}}^{K,{(k_{j})}}} \right)}\mspace{14mu}\ldots}} \\{p\left( {x_{n}^{P_{j},{(i_{j})}}❘x_{n - 1}^{P_{j},{(k_{j})}}} \right)} \\{\approx {\underset{k_{j}}{argmax}\alpha_{j,{n - 1}}^{(k_{j})}{p\left( {{y_{n}❘x_{n}^{P_{j},{(i_{j})}}},{\hat{x}}_{n❘{0:{n - 1}}}^{K,{(i_{j})}}} \right)}\mspace{14mu}\ldots}} \\{p\left( {x_{n}^{P_{j},{(i_{j})}}❘x_{n - 1}^{P_{j},{(k_{j})}}} \right)} \\{= {\underset{k_{j}}{argmax}\alpha_{j,{n - 1}}^{(k_{j})}{{p\left( {x_{n}^{P_{j},{(i_{j})}}❘x_{n - 1}^{P_{j},{(k_{j})}}} \right)}.}}}\end{matrix} & \begin{matrix}(42) \\\; \\\; \\\; \\(43) \\\; \\\; \\\; \\(44)\end{matrix}\end{matrix}$

Given the best trajectory for each particle i_(j), corresponding Kalmanstate variables χ^(P) ^(j,) ^((i) ^(j) ⁾, i.e. sinusoidal coefficients,are updated utilizing the extended Kalman filter. Then, the MAP estimateof χ_(n) ^(P) ^(i) is obtained based on the value of the coefficientα_(j,n) ^((i) ^(j) ⁾ as follows,

$\begin{matrix}{i_{j,n}^{*} = {\arg{\max\limits_{i_{j}}\alpha_{j,n}^{(i_{j})}}}} & (45) \\{{\hat{x}}_{n}^{P_{i}} = x_{n}^{P_{i},{(i_{j,n}^{*})}}} & (46)\end{matrix}$

Since there are two groups of particles i₁ and i₂, we need to select thebest estimate of the Kalman state vector χ_(n) ^(K) among two potentialestimates: χ_(n) ^(K,(i) ^(1,n) ^(*)) and χ_(n) ^(K,(i) ^(2,n) ^(*)).The actual estimate of the Kalman state vector χ_(n) ^(K) can beselected as follows,

$\begin{matrix}{i_{{MAP},n}^{*} = \left\{ {\begin{matrix}i_{1,n}^{*} & {\alpha_{1,n}^{(i_{1,n}^{*})} \geq \alpha_{2,n}^{(i_{2,n}^{*})}} \\i_{2,n}^{*} & {\alpha_{1,n}^{(i_{1,n}^{*})} < \alpha_{2,n}^{(i_{2,n}^{*})}}\end{matrix}.} \right.} & (47) \\{{\hat{x}}_{n}^{K} = x_{n}^{K,{(i_{{MAP},n}^{*})}}} & (48)\end{matrix}$

Then, the estimate of the entire state χ_(n) can be expressed as,{circumflex over (χ)}_(n)={{circumflex over (χ)}_(n) ^(P) ^(1,) ^((i)^(1,n) ^(*)),{circumflex over (χ)}_(n) ^(P) ^(2,) ^((i) ^(2,n)^(*)),{circumflex over (χ)}_(n) ^(K,(i) ^(MAP,n) ^(*))}.  (49)

Algorithm 2 Table - Dual MAM-PF Method Steps Initializationfor  j = 1, 2  do   for  i = 1, …  , N_(p)  do    Sample  x₀^(P_(j), (i)) ∼ π₀(x₀^(P_(j)))    x̂_(0 : −1)^(K, (i)) = E[x₀^(K, (i))|x₀^(P, (i))]    α₀^((i)) = π₀(x₀^(P, (i)))p(y₀|x₀^(P, (i)), x₀^(K, (i)))    z₀^((i)) = x₀^((i))   end  for$\mspace{25mu}{i_{j,0}^{*} = {\underset{i_{j}}{argmax}\mspace{14mu}\alpha_{j,o}^{(i_{j})}}}$end  for$i_{{MAP},0}^{*} = {\underset{i_{j,0}^{*}}{argmax}\mspace{14mu}\alpha_{j,0}^{(i_{j,0}^{*})}}$x̂₀ = {x₀^(P_(1,)(i_(1, 0)^(*))), x₀^(P_(2,)(i_(2, 0)^(*))), x₀^(K_(,)(i_(MAP, 0)^(*)))}for  n = 1, …  , N_(T)  do   for  j = 1, 2  do    for  i_(j) = 1, …  , N_(p)  do      MAP  Estimation       x_(n)^(P_(j), (i_(j))) ∼ q_(n)(x_(n)^(P_(j), (i_(j)))|x_(n − 1)^(P_(j), (i_(j))), y_(n))$\mspace{101mu}{k_{j}^{*} = {\underset{k_{j}}{argmax}\mspace{14mu}\alpha_{j,{n - 1}}^{(k_{j})}{p\left( x_{n - 1}^{P_{j},{(i_{j})}} \middle| x_{n - 1}^{P_{j},{(k_{j})}} \right)}}}$       Marginalized  Sequential  Estimation       x_(n)^(P, (i_(j))) = {x_(n)^(P₁, (i_(1, n − 1)^(*))), …  , x_(n)^(P_(j), (i_(j))), …  , x_(n)^(P_(N_(j), (i_(N_(j), n − 1)^(*))))}       ŷ_(n|0 : n − 1) = H_(n)(x_(n)^(P, (i_(j))))x̂_(n|0 : n − 1)^(K, (k_(j)^(*)))       e_(n) = y_(n) − ŷ_(n|0 : n − 1)       R_(v, n) = [e_(n)e_(n)^(T) − H_(n)(x_(n)^(P, (i_(j))))C_(n|0 : n − 1)^((k_(j)^(*)))H_(n)(x_(n)^(P, (i_(j))))^(T)]₊      R̂_(v, n) = β R̂_(v, n − 1)^((k_(j)^(*))) + (1 − β)R_(v, n)      R_(e, n) = H_(n)(x_(n)^(P, (i_(j))))C_(n|0 : n − 1)^((k_(j)^(*)))H_(n)(x_(n)^(P, (i_(j))))^(T) + R̂_(v, n)       K_(n) = C_(n|0 : n − 1)^((k_(j)^(*)))H_(n)(x_(n)^(P, (i_(j))))^(T)(R_(e, n))⁻¹       x̂_(n|0 : n)^(K, i_(j)) = x̂_(n|0 : n − 1)^(K, (k_(j)^(*))) + K_(n)e_(n)       C_(n|0 : n) = [I − K_(n)H_(n)(x_(n)^(P, (i_(j))))]C_(n|0 : n − 1)^((k_(j)^(*)))       x̂_(n + 1|0 : n)^(K, (i_(j))) = F_(n)(x_(n)^(P, (i_(j))))x̂_(n|0 : n)^(K, i_(j))&  R̂_(v, n)^((i_(j))) = R̂_(v, n)       C_(n + 1|0 : n)^((i_(j))) = F_(n)(x_(n)^(P, (i_(j))))C_(n|0 : n)F_(n)(x_(n)^(P, (i_(j))))^(T) + Q_(u)^(K)       α_(j, n)^((i_(j))) = α_(j, n − 1)^((k_(j)^(*)))p(y_(n)|x_(n)^(P, (i_(j))), x̂_(n|0 : n − 1)^(K, (k_(j)^(*))))        p(x_(n)^(P_(j), (i_(j)))|x_(n − 1)^(P_(j), (k_(j)^(*))))     end  for      Update  MAP  State  Estimate$\mspace{79mu}{i_{j,n}^{*} = {\underset{i_{j}}{argmax}\alpha_{j,n}^{(i_{j})}}}$    end  for$\mspace{50mu}{i_{{MAP},n}^{*} = {\underset{i_{j,n}^{*}}{argmax}\mspace{14mu}\alpha_{j,n}^{(i_{j,n}^{*})}}}$    x̂_(n) = {x̂_(n)^(P₁, (i_(1, n − 1)^(*))), x̂_(n)^(P₂, (i_(2, n)^(*))), x̂_(n)^(K, (i_(MAP, n)^(*)))}end  forH.4. ABP Signal Envelope Estimation

Given the estimated signal parameters in (26)-(34), according to oneembodiment, it is possible to estimate the upper envelope (e_(μ,n)) andlower envelope (e_(l,n)) of ABP signals by following steps below,

$\begin{matrix}{{\theta_{\max,n}^{c} = {\arg\;{\max\limits_{\theta}{\sum\limits_{k = 1}^{N_{h}^{c}}{\rho_{k,n}\left\lbrack {{c_{1,k,n}{\cos\left( {k\;\theta} \right)}} + {c_{2,k,n}{\sin\left( {k\;\theta} \right)}}} \right\rbrack}}}}}{\theta_{\min,n}^{c} = {\arg\;{\min\limits_{\theta}{\sum\limits_{k = 1}^{N_{h}^{c}}{\rho_{k,n}\left\lbrack {{c_{1,k,n}{\cos\left( {k\;\theta} \right)}} + {c_{2,k,n}{\sin\left( {k\;\theta} \right)}}} \right\rbrack}}}}}{\kappa_{\max,k,n} = {{c_{1,k,n}{\cos\left( {k\;\theta_{\max,n}^{c}} \right)}} + {c_{2,k,n}{\sin\left( {k\;\theta_{\max,n}^{c}} \right)}}}}{\kappa_{\min,k,n} = {{c_{1,k,n}{\cos\left( {k\;\theta_{\min,n}^{c}} \right)}} + {c_{2,k,n}{\sin\left( {k\;\theta_{\min,n}^{c}} \right)}}}}} & \; \\{e_{\mu,n} = {\gamma_{n} + {\sum\limits_{k = 1}^{N_{h}^{c}}{\rho_{k,n}\kappa_{\max,k,n}}}}} & (50) \\{e_{l,n} = {\gamma_{n} + {\sum\limits_{k = 1}^{N_{h}^{c}}{\rho_{k,n}\kappa_{\min,k,n}}}}} & (51)\end{matrix}$where arg max_(χ)ƒ(χ) and arg min_(χ)ƒ(χ) are operators to obtain thevalue of χ for which ƒ(χ) attains its maximum and minimum values,respectively.H.5. Pulse Pressure Signal Envelope Estimation

The pulse pressure (PP) signal is the difference between the upperenvelope e_(μ,n) and lower envelope e_(l,n) of the ABP signal. This PPsignal oscillates roughly at the respiratory rate. This oscillation isdue to the respiratory effect on the variation of systemic ABP underfull ventilation support. Within each respiratory cycle PP reaches itsmaximum (PP_(max)) and minimum (PP_(min)) values, which are two criticalparameters to compute the PPV index. Traditionally, the PP_(max) andPP_(min) values have been computed only once per each respiratory cycle.Given the estimated signal parameters in (26)-(34), however, thecontinuous equivalents of PP_(max) and PP_(min) can be computed. Theyare the upper ε_(μ,n), and lower ε_(l,n) envelopes of the PP signal. Theupper envelope ε_(μ,n) is the continuous estimate of PP_(max) and thelower envelope ε_(l,n) that of PP_(min). The ε_(μ,n) and ε_(l,n) valuescan be estimated as described below,

$\varrho_{k,n} = {{\sum\limits_{j = 1}^{N_{h}^{r}}{m_{1,k,j,n}{\cos\left( {j\;\theta} \right)}}} + {m_{2,k,j,n}{\sin\left( {j\;\theta} \right)}}}$$\theta_{\max,n}^{r} = {\arg\;{\max\limits_{\theta}{\sum\limits_{k = 1}^{N_{h}^{c}}{\left( {1 + \varrho_{k,n}} \right)\left( {\kappa_{\max,k,n} - \kappa_{\min,k,n}} \right)}}}}$$\theta_{\min,n}^{r} = {\arg\;{\min\limits_{\theta}{\sum\limits_{k = 1}^{N_{h}^{c}}{\left( {1 + \varrho_{k,n}} \right)\left( {\kappa_{\max,k,n} - \kappa_{\min,k,n}} \right)}}}}$$\varrho_{\max,k,n} = {{\sum\limits_{j = 1}^{N_{h}^{r}}{m_{1,k,j,n}{\cos\left( {j\;\theta_{\max,n}^{r}} \right)}}} + {m_{2,k,j,n}{\sin\left( {j\;\theta_{\max,n}^{r}} \right)}}}$$\varrho_{\min,k,n} = {{\sum\limits_{j = 1}^{N_{h}^{r}}{m_{1,k,j,n}{\cos\left( {j\;\theta_{\min,n}^{r}} \right)}}} + {m_{2,k,j,n}{\sin\left( {j\;\theta_{\min,n}^{r}} \right)}}}$$\begin{matrix}{ɛ_{\mu,n} = {\sum\limits_{k = 1}^{N_{h}^{c}}{\left( {1 + \varrho_{\max,k,n}} \right)\left( {\kappa_{\max,k,n} - \kappa_{\min,k,n}} \right)}}} & (52) \\{ɛ_{l,n} = {\sum\limits_{k = 1}^{N_{h}^{c}}{\left( {1 + \varrho_{\min,k,n}} \right)\left( {\kappa_{\max,k,n} - \kappa_{\min,k,n}} \right)}}} & (53)\end{matrix}$where 1+

is equal to ρ_(k,n), and ε_(μ,n) and ε_(l,n) are the continuousestimates of the PP_(max) and PP_(min), respectively.H.6. Pulse Pressure Variation Calculation

Given the ε_(μ,n) and lower ε_(l,n) values, it is possible to calculatethe PPV index. According to one embodiment, it can be computed asfollows,

$\begin{matrix}{{P\; P\;{V(\%)}} = {100 \times \frac{ɛ_{\max} - ɛ_{\min}}{\left( {ɛ_{\max} + ɛ_{\min}} \right)/2}}} & (54)\end{matrix}$

This new PPV index is different from the traditional PPV index describedin (1) because the new one is continuous in time while the traditionalone can be obtained only once per each respiratory cycle.

While particular embodiments have been described, it is understood that,after learning the teachings contained in this disclosure, modificationsand generalizations will be apparent to those skilled in the art withoutdeparting from the spirit of the disclosed embodiments. It is noted thatthe foregoing embodiments and examples have been provided merely for thepurpose of explanation and are in no way to be construed as limiting.While the apparatus has been described with reference to variousembodiments, it is understood that the words which have been used hereinare words of description and illustration, rather than words oflimitation. Further, although the system has been described herein withreference to particular means, materials and embodiments, the actualembodiments are not intended to be limited to the particulars disclosedherein; rather, the system extends to all functionally equivalentstructures, methods and uses, such as are within the scope of theappended claims. Those skilled in the art, having the benefit of theteachings of this specification, may effect numerous modificationsthereto and changes may be made without departing from the scope andspirit of the disclosed embodiments in its aspects.

The invention claimed is:
 1. A method for assessing fluid responsivenesscomprising: (a) obtaining a cardiovascular signal; and (b) computing adynamic index predictive of fluid responsiveness from saidcardiovascular signal using a hardware device with one or moreprocessors configured for implementing a nonlinear state spaceestimator, said nonlinear state space estimator comprising a state spacemodel wherein a plurality of cardiac harmonics are separately modulatedby respiratory fluctuations; whereby said dynamic index predictive offluid responsiveness enables a practitioner to assess fluidresponsiveness and guide fluid therapy.
 2. The method of claim 1,wherein said cardiovascular signal is an arterial blood pressure signal.3. The method of claim 1, wherein said cardiovascular signal is aplethysmographic signal.
 4. The method of claim 1, wherein said dynamicindex predictive of fluid responsiveness is a nonlinear function of saidcardiovascular signal.
 5. The method of claim 4, wherein said dynamicindex predictive of fluid responsiveness is substantially equivalent toa variation in pulse pressure of said cardiovascular signal due torespiration.
 6. The method of claim 5, wherein said nonlinear statespace estimator is based on a nonlinear state space model ofcardiovascular signals.
 7. The method of claim 6, wherein said nonlinearstate space estimator is a computational Bayesian estimator.
 8. Themethod of claim 7, wherein said computational Bayesian estimator is asequential Monte Carlo estimator.
 9. The method of claim 8, wherein saidsequential Monte Carlo estimator is a marginalized particle filter. 10.The method of claim 9, wherein said marginalized particle filter is ana-posteriori adaptive marginalized particle filter.
 11. An apparatus forassessing fluid responsiveness, comprising: (a) a measuring unit formeasuring a cardiovascular signal; and (b) a processor configured forcomputing a dynamic index predictive of fluid responsiveness from saidcardiovascular signal using a nonlinear state space estimator, saidstate space estimator comprising a state space model wherein a pluralityof cardiac harmonics are separately modulated by respiratoryfluctuations; whereby said dynamic index predictive of fluidresponsiveness enables a practitioner to assess fluid responsiveness andguide fluid therapy.
 12. The apparatus of claim 11, wherein said dynamicindex predictive of fluid responsiveness is substantially equivalent toa variation in pulse pressure of said cardiovascular signal due torespiration.
 13. The apparatus of claim 12, wherein said nonlinear statespace estimator is based on a nonlinear state space model of saidcardiovascular signal.
 14. The apparatus of claim 13, wherein saidnonlinear state space estimator is a computational Bayesian estimator.15. The apparatus of claim 14, wherein said computational Bayesianestimator is a sequential Monte Carlo estimator.
 16. The apparatus ofclaim 15, wherein said sequential Monte Carlo estimator is amarginalized particle filter.
 17. The apparatus of claim 16, whereinsaid marginalized particle filter is an a-posteriori adaptivemarginalized particle filter.
 18. The apparatus of claim 17, whereinsaid marginalized particle filter is adapted to work on a cardiovascularsignal from a spontaneously breathing patient.